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In this paper, we survey five different computational modeling methods. For comparison, we use 
the activation cycle of G-proteins that regulate cellular signaling events downstream of G-protein- 
coupled receptors (GPCRs) as a driving example. Starting from an existing Ordinary Differential 
Equations (ODEs) model, we implement the G-protein cycle in the stochastic Pi-calculus using 
SPiM, as Petri-nets using Cell Illustrator, in the Kappa Language using Cellucidate, and in Bio- 
PEPA using the Bio-PEPA eclipse plug in. We also provide a high-level notation to abstract away 
from communication primitives that may be unfamiliar to the average biologist, and we show how to 
translate high-level programs into stochastic Pi-calculus processes and chemical reactions. 



1 Introduction 



Traditionally, biologists have used ordinary differential equations (ODEs) to model biological processes 
and simulate the evolution of species over time. However, the past decade has seen the emergence of a 
family of formalisms dedicated to modeling biological behavior (kinetics). Stemming from formal lan- 
guages designed to capture concurrent computation and communicating processes, these new formalisms 
range from graphical formalisms to algebraic languages. 

In this paper we survey five different computational modeling formalisms, and we show how to 
simulate the activation cycle of G-proteins by G-protein coupled receptors in each of them. 



1.1 The G-protein cycle 

G-protein-coupled receptors (GPCRs) are seven-pass transmembrane receptors ll23l that mediate extra- 
cellular signaling molecules and intracellular signal transduction pathways via trimeric GTP binding 
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proteins (G proteins) ifTTIl . It is estimated that over half of all marketed pharmaceuticals target GPCRs 
EH. 

G-protein-coupled receptor's primary function is to transduce extracellular stimuli into intracellular 
signals El[T6]|(Fig.[T]). These receptors' intracellular domains interact with heterotrimeric G-proteins Q. 
G-protein-coupled receptors are commanded both by a wealth of stimuli to which they respond, as well 
as by the assortment of intracellular signaling pathways they activate. These include light, neurotrans- 
mitters, odorants, biogenic amines, lipids, ligands, hormones, nucleotides, and chemokines fl6l . 

G-proteins are composed of a, j6 , and y subunits. ft and y form the j8 7 subcomplex, and the a subunit 
can be bound to GDP or GTP. When a ligand activates the G-protein-coupled receptor, it induces a 
conformational change in the receptor that exchanges GDP for GTP on the a subunit and this triggers 
the dissociation of the a subunit, which is bound to GTP, from the /3y dimer and the receptor El[T6l. 
Then the a subunit will act on the target proteins. After that, the bound GTP will be hydrolyzed to GDP, 
and the a unit, which is bound to GDP, will bind the j8 7 dimer and the receptor again. From step 5 to step 
l(Fig. [TJ, hydrolysis from GTP to GDP happens, and is enhanced by binding a regulator of G-protein 
signaling(RGS) MM- 




Fi gure 1 '. Activation cycle of G-proteins by G-protein-coupled receptors. When a ligand activates the receptor, a conformation 
change occurs in the receptor that exchanges GDP for GTP on the a subunit and this triggers the dissociation of the a subunit 
from P 7 dimer and the receptor. After the free a unit works on target proteins, GTP will be hydrolyzed to GDP. The GTPase 
activity is enhanced by binding of the RGS. 

While the mechanism of G-protein signaling is highly conserved, the diversity of signaling targets is 
immense | 26 |. Thus, the ability to swiftly model the dynamics of G-protein signaling is advantageous to 
both drug development and basic research. 

1.2 Modeling 

Two recent studies use ODEs to model the receptor and ligand binding dynamics and the G-protein 
cycle. Yi. et. al. |[28l used an ODEs model together with fluorescent resonance energy transfer (FRET) 
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measurement to characterize the heterotrimeric G-protein cycle in yeast. Hao, et. al. lfT2l also used an 
ODEs model to indicate the role of Sst2, an RGS protein, in G-protein cycling. 

In our study, we build four new models using the stochastic Pi-calculus, stochastic Petri Nets, the 
Kappa language, and Bio-PEPA. Furthermore, we show that we obtain consistent results when compared 
with the ODEs model of Yi, et. al. lf28l (Fig. [5]). The simulators for stochastic Petri Nets ll20l and the 
Kappa language iflOlL Cell Illustrator and Cellucidate, respectively, have an interface which is intuitive for 
biologists, and the syntax of Bio-PEPA is close to the chemical reactions notation and relatively easy to 
use, unlike SPiM, which requires understanding of communications primitives and concurrent processes. 
In order to hide such communication details we develop a high level notation that uses terminology 
directly obtained from biological processes, and we show how to systematically translate a process in 
high level notation into a SPiM process. (See Section |4j). Kahramanogullari. et. al. llT5l and Guerriero. 
et. al. [llj have provided the narrative languages for SPiM and Bio-PEPA, respectively. Our high level 
notation is an alternative to their narrative languages. 

1.3 Outline 

In Section [2j we show how to implement the model for the G-protein cycle using five different meth- 
ods. In Section [3] we compare the simulation results from the different implementations. In Section |4| 
we introduce a high level notations and the translations from high level notation to chemical reactions 
and SPiM code. Section [5] contains our conclusions and a comparison of the simulation and modeling 
approaches of the five methods considered in Section [2] 

2 Method 

2.1 ODEs modeling 

The ODEs modeling approach is a traditional method applied by biologists to explore dynamic biological 
processes by numerical integration. Yi. et. al. ll28l give a computational model of the G-Protein cycle by 
means of ordinary differential equations (ODEs). The individual reactions comprising the key dynamics 
of heterotrimeric G-proteins in yeast are represented, along with their rate constants, in Table [T] 

Table 1 : The reactions used to model the G-protein cycle and the corresponding rate constants (rate parameters) for each 
reaction (G-protein tutorial from www.MathWorks.com). 



No. 


Reaction 


Rate Parameters 


1 


L + R o RL 


Krl = 2- \O b M~ l -s-\ K RLm = 1 • 10" V 1 


2 


Gd + Gbg ->■ G 


Kgi = Is -1 


3 


RL + G -»■ Ga + Gbg + RL 


K Ga = \- 10-V 1 


4 


R o null 


K Rd o = 4-10- 4 S - l ,K Rs = 4s- 1 


5 


RL ->• null 


^^4-10-V 1 


6 


Ga-> Gd 


K G(i = 0.11s" 1 



Reaction No. 1 is receptor-ligand interaction which corresponds to the transition from step 2 to step 
3 in Fig. [T] R represents receptor; L represents ligand; and RL represents receptor and ligand bound 
together. Reaction No. 2 is heterotrimeric G-protein formation which can be found in the transition 
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from step 1 to step 2 in Fig. [T] Gd represents the deactivated a subunit, which is bound to GDP; Gbg 
represents the /3 y sub-complex, and G represents a bound to /3y and GDP. Notice that from step 1 to 2, 
the a subunit also binds to the receptor. However, the original paper does not consider it as a species, 
and consequently the ODEs model does not reflect that binding. Reaction No. 3 is G-protein activation 
which corresponds to the transition from step 3 to step 4 in Fig. [T] After the ligand and receptor interact 
with each other, the a subunit will be activated and bound to GTP, and it will dissociate from the /3y 
sub-complex. Ga represents the activated a subunit, which is bound to GTP. Reaction No. 4 is receptor 
synthesis and degradation. Reaction No. 5 is receptor-ligand degradation. Reaction No. 6 is G-protein 
inactivation which corresponds to the transition from step 5 to step 1 in Fig.[T| 
The ODEs of the model are as follows: 

i. d[R]/dt=-K RL [L][R] + K RLm [RL]-K Rd0 [R] +K Rs ; 

ii. d[RL]/dt=K RL [L] [R] - K RLm [RL] - K Rdl [RL] ; 
hi. d[G]/dt=-K Ga [RL] [G] + K G1 [Gd] [Gbg] ; 

iv. d[Ga]/dt=K Ga [RL] [G] - K Gd [Ga] ; 

v. d[Gd]/dt=-K Gi [Gd] [Gbg] + K Gd [Ga] ; 

vi. d[Gbg]/dt=K Ga [RL] [G] - K G1 [Gd] [Gbg] ; 

vii. d[L]/dt=-K RL [L] [R] + K RLm [RL] ; 

According to Table[T] derived from the reactions diagram from Yi. et. al. (281, there are seven species 
in the model. However, they only describe the first four equations. In order to have one differential equa- 
tion for each species, we build the last three equations. They also provide two conservation relationships: 
[Gbg] = Gt — [G] and [Gd] = Gt — [G] — [Ga] which are satisfied by our simulations. (See Section[3]). 

2.2 Stochastic Pi-calculus modeling 

Process algebras are formal languages originally designed to model complex reactive computer sys- 
tems where heterogeneous agents interact concurrently exchanging information through communication 
channels. Because of the similarities between reactive computer systems and biological systems, process 
algebras have recently been used to model biological systems Il6ll25ll. 

The stochastic Pi-calculus 125] is a process algebra where stochastic rates are imposed on processes, 
allowing for more accurate description of biological systems |6|. A process can be depicted as a collec- 
tion of interacting automata with two kinds of reactions: delay @r and interaction® r on ch. A delay is a 
spontaneous change of state performed by a single automaton at a specified reaction rate r. Interaction, 
on the other hand, consists of a send/receive handshake between two automata over a channel ch, written 
!ch for send and ?ch for receiveQ The channel ch has an associated reaction rate r. Reaction rates are 
used to determine stochastically the reaction to be executed next [5|. For our modeling, we use SPiM, an 
implementation of the stochastic Pi-calculus that can be used to run in-silico simulations displaying the 
change in the populations of different species over time ll22l . 

Fig. [2] depicts the stochastic Pi-calculus model we build for the G-protein cycle using SPiM. We 
divide the system into two parts that correspond to two viewpoints. The left hand side is the viewpoint 
of ligand and receptor, and the right hand side is the viewpoint of the G-protein. The initial numbers of 
different species are assigned in square brackets. For example, the initial concentration of R is 10000 and 



Note that it is sufficient to use the simplest form of communication where !ch and ?ch exchange an empty message. 
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[10000] 



(^lT^) [602200] 




Figure 2: Graphical representation of Pi-calculus modeling 



of G is 7000. Notice that RL and Ga are not initialized. All the reaction rates are also extracted from the 
literature ll28l and continuous reactions are converted to discrete reactions according to standard rules J?). 
For example, the pair !bindR/?bindR represents step 2 to 3 in Fig.[T} and it corresponds to reaction No. 1 
in Table [TJ and the pair !switch/?switch represents the step 3 to 4 in Fig.[T] and it corresponds to reaction 
No. 3 of Table[T} The corresponding SPiM code is in Table|2] The "directive plot" statement declares the 
species that will be plotted, in this case RL, G, and Ga. At the end of the program, each "run" statement 
declares the initial amount for different species. In the ODEs model, the amount of R is initialized to 
be 6.022E17 and the rate of the reaction No. 1 in Table[T]is 3.32E-18 with the assumed reaction volume 
1 L (M = mol L _1 , mol = 6.022 x 10 23 ). Because simulating stochastic models is slower than solving 
ODEs, we scale the reaction volume down to 1.0E-12 L in order to speed up the stochastic simulation. 
Given this new reaction volume, the initial amount of R is 6.022E5, and the new reaction rate is 3.32E- 
6|@1. The other reaction rates remain unchanged, because they are discrete and do not depend on the 
volume. Assuming that pressure and temperature are unchanged, proportionally scaling reaction rates 
and reagent concentrations preserves the correctness of the modeling. We do not model the receptor 
generation aspect of reaction No. 4 in Table [T] (null —> R), because on the one hand, this reaction is 
inconsequential to the whole system, as suggested by the result of Fig. |5| and on the other hand, the 
stochastic pi-calculus can not describe the generation of a process from the empty process. Furthermore, 
the initial amounts of receptor and G-protein guarantee abundance of receptor to enable the activation of 
G-protein. More details about our model can be found in the companion technical report Q. 



2.3 Petri Nets modeling 

A Petri Net is a mathematical concept developed in the early 1960s by Carl Adam Petri to describe 
discrete distributed systems. Petri Nets are suitable for modeling and analyzing the dynamics of con- 
current systems whose behavior could be described by finite sets of atomic processes and atomic states 
ll20l . Petri Nets have subsequently been adapted and extended in many directions, including quantitative 
analysis of biological networks ll24l . 

The basic Petri Net is a directed bipartite graph with two kinds of nodes which are either places 
or transitions and directed arcs which connect nodes. In modeling biological processes, place nodes 
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Table 2: SPiM code for G-protein Cycle 



directive sample 600.0 
directive plot RL(); G(); Ga() 

val rl = 3.32E-6 
val rlm= 0.01 
val rs= 4.0 
val rd0= 4.0E-4 
val rdl= 4.0E-3 
val ga= 1.0E-5 
val gdl= 0.11 
valgl= 1.0 

new bindb@ 



new bindR@rl:chan 
new switch @ga:chan 

let R()= 

( 

do ?bindR(); RL() 
or delay @rd0; () 

) 

and RL()= 

( 

do delay @ rim; ( R() 
or delay @rdl; () 
or !switch(); RL() ) 



L0) 



and L()=( !bindR(); () ) 

let Gd()= ( !bindb(); G() ) 
and GO ( ?switch(); (Ga() | Gbg() ) ) 
and Ga()= ( delay @gdl; Gd()) 
and Gbg()= ( ?bindb(); () ) 

run 602200 of L() 
run 10000 of R() 
run 7000 of G() 
run 3000 of Gd() 
run 3000 of Gbg() 



represent molecular species and transition nodes represent reactions. We use Cell Illustrator to develop 
our model of the G-protein cycle ( www . cellillustrator . co"m| ). 

Fig. [3] is the Petri Net we build to model the G-protein cycle. Place nodes are represented as circles, 
and transition nodes are represented as boxes. The arcs are labeled with an integer weight, which rep- 
resents the minimum value of input entities needed for the reactions to happen. The initial number of 
species and reaction rates can be read from the graph. For example, the initial amount of Receptor is 
10000 and the initial amount of j8 y sub-complex is 3000, while Ga and RL are initialized to 0. Here we 
scale the reaction volume to be 1.0E-12 L, which is the same as in the SPiM model. The numbers in red 
correspond to the No. column in Table [T] 

2.4 Kappa language modeling 



Table 3: G-protein cycle is represented in Kappa language. 



No. 


Kappa Statement 


1 


R(r), L(l) o R(r!l), L(l!l) @ 3.32e-06,0.01 


2 


Gbg(bg), Gd(alpha - GDP,a) Gbg(bg!l), Gd(alpha - GDP,a!l) @ 1.0 


3 


Gbg(bg!l), Gd(alpha - GDP,a!l), R(r!2), L(l!2) Gbg(bg), 
Gd(alpha - GTP,a), R(r!l), L(l!l) @ 1.0e-05 


4 


R(r)++ @4e-04,4 


5 


L(l!l), R(r!l)-> @ 0.004 


6 


Gd(alpha - GTP,a) Gd(alpha - GDP,a) @ 0.11 



Kappa is a formal language for defining agents (typically meant to represent proteins) as sets of 
sites that constitute abstract resources for interaction. It is used to express rules of interactions between 
proteins characterized by discrete modification and binding states fLOl . The Kappa language modeling 
platform Cellucidate is available at www . cellucidate . com) 

In the Kappa language, reaction rules are described by rewriting rules between lists of agents. Each 
agent has a name and binding sites. For example, in R(r), R is the agent and r is its binding site. Agents 
can become bound, and the two end points of a link between two agents is indicated by !i, for some 
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Figure 3: Petri Nets modeling 

index i. rvalue specifies the internal state of a site on the agent, <H> specifies a bidirectional reaction, —> 
specifies an unidirectional reaction, and ©value specifies the reaction rate. For example, 

R(r), L(l) R(r!l), L(l!l) @ 3.32e-06,0.01 

is a bidirectional reaction, where R and L are agent names, r is the binding site of R, 1 is the binding 
site of L; 0.01 is the reaction rate from left to right, and 3.32e-06 is the reaction rate from right to left. 
In R(r! 1) and L(l! 1), 1 is the index of the link that binds R and L at their binding sites. In our G-protein 
example, it corresponds to the binding of receptor and ligand. Gd( a^GDP,a) represents the deactivated 
G-protein. ^GDP represents the internal state of the a subunit where it is bound to GDR 

The Kappa language processes for the G-protein cycle are described in Table [3] The statements 
correspond to the reactions and rates in Table The initial number of species can not be read from 
Table |3j but in fact, we use the same number as the Petri Nets model and SPiM model. For example 
the initial number for ligand is 6.022E5. Fig.|4]is the graphical representation of the processes described 
in Table [3] automatically generated by Cellucidate. R + LRL means R + L —> RL, and R + LRL_op is 
the opposite direction: R + L ^— RL. NullR means R —¥ null, and NullR^op is the opposite direction: 
R ^— null. The reactions with red arrows will consume R or RL, having a negative effect on the reaction 
R + RL^ Ga + Gbg + RL, while the green arrows have positive effects. 

2.5 Bio-PEPA modeling 

Bio-PEPA[9| is an extension of PEPA lfT3lK which is specifically designed for modeling biochemical net- 
works by explicitly defining the stoichiometry of reactions. The two existing tools for working with 
Bio-PEPA models are the Bio-PEPA Workbench and the Bio-PEPA Eclipse Plug-in||8]|. We use the Bio- 
PEPA Eclipse Plug-in to model the G-protein Cycle, shown in Table |4j In Table |4| Parameter Definitions 

2 Note that the initial amount of different species can be set in Cellucidate. 
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Fi gure 4: Graphical representation of Kappa language modeling 

consist of rate constant declarations; Functional Rates are declarations of reaction names and the corre- 
sponding reaction rates. For example, in kineticLowOfbindRL : fMA(rl), the name of the reaction is 
bindRL, the function for the law of Mass-Action is denoted by fMA, and rl is the rate constant; under 
Species Components there is a representation of the reactions matching the corresponding differential 
equations of the ODEs model. For example, 

G = bindAB » G + activ « G 

indicates the role of reagent G in the related reactions. bindAB » G means G is produced in reaction 
bindAB and activ «G means G is consumed in reaction activ. It also corresponds to 

d[G]/dt=-K Ga [RL] [G] + Kg i [Gd] [Gbg] , 
where gl is the rate constant of reaction bindAB, and ga is the rate constant of reaction activ, corre- 
sponding to Kgi and K Ga respectively. 

Finally, the initial concentrations for each species are declared under Model Components. Because 
in our example the reactions are in the same cell, and the locations are the same for all species, the default 
location is used and left implicit. 

The reaction rates and initial concentrations are all consistent with the other three models we built. 

3 Simulation Results and Comparison 

Although we are not proposing a theoretical comparison of the five modeling techniques, or the five im- 
plementations in different formalisms, we show here consistent experimental simulation results. SPiM, 
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Table 4: Bio-PEPA code for G-protein Cycle 

// Parameter Definitions 

r/ = 3.32£-6; 

rlm = 0.01; 

ga= 1.0£-5; 

gdl =0.11; 

gl = 1.0; 

rd\ = 4.0£-3; 

= 4.0; 
rd0 = 4.0£-4; 

// Functional Rates 
kineticLawOfbindRL : fMA(rl)\ 
kineticLawOfbindAB : /MA(gl); 
kineticLawOfdisRL : fMA(rlm)\ 
kineticLawOfactiv : fMA(ga); 
kinetic Law Of hydro : fMA(gdl); 
kineticLawOfdegRL : fMA(rd\)\ 
kineticLawOfdegR : /MA(rdO); 
kineticLawOfproR : as; 

// Species Components 

L = dwftL > > L + bindRL «L\ 

R = disRL »R + bindRL «R + degR «R + proR » R\ 

RL = disRL «RL + bindRL » RL + activ(+)RL + degRL < < RL\ 

G = bindAB » G + activ « G; 

Ga = art/v >> G<2 + hydro << Ga\ 

Gd bindAB «Gd + hydro » Gd\ 

Gbg bindAB « Gbg + activ » Gbg\ 

II Model Components 

L[602200] < * > #[10000] < * > RL[0] < * > G[7000] 

< * > G d [3000] < * > Gbg[3000] < * > Ga[0] 



Cell Illustrator, Cellucidate, and Bio-PEPA all use Gillespie's algorithm for simulation. Fig. [5] shows 
that the results from those four simulations are consistent with the result of the original ODEs modeling. 
The plots show that RL is consumed as the reaction proceeds. The curves of G and Ga decline and 
increase oppositely, and the sum of these two species is a constant, which is equal to the initial amount 
of G-protein, so the conservation relationships mentioned by Yi. et. al. are confirmed ([Gbg] =Gt— [G], 
[Gd] = Gt-[G]-[Ga])tiM. 

These computational models can be classified according to different properties. 

Simulation Approach. The simulation approach can be stochastic or deterministic. While the de- 
terministic approach is faster, the stochastic approach is more accurate. In particular, because of the 
randomness of dynamic biological systems, the stochastic approach can make the simulation of biolog- 
ical processes more faithful. Furthermore, while a deterministic system may get stuck into a specific 
state, stochastic noise allows the system to oscillate in and out of that state BTll . 

Our four models of the G-protein cycle in the stochastic Pi-calculus, stochastic Petri Nets, the Kappa 
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language, and Bio-PEPA fall into the category of stochastic simulation, while the traditional ODEs model 
is deterministic. 

Far from competing with each other, both modeling approaches are complementary, and often it is 
possible to derive one model from the other. In particular, Feret, et. al. show how to convert a rule-based 
model, such as a Kappa model, into a reduced systems of differential equations fLOl . Priami, et. al. show 
how to translate an ODEs model into a Stochastic Pi-calculus model fl9lk and Cardelli, et. al. show how 
to convert processes to ODEs 01. 

Modeling Approach. The five models discussed here are either chemical reaction centric or process 
centric. The ODEs, Petri Nets, Kappa, and Bio-PEPA models are chemical reaction centric. In particular, 
Kappa takes a rule based approach that corresponds to high level chemical reactions. For each chemical 
reaction there is an equation, a transition or a rule in the corresponding model. On the other hand, the 
stochastic Pi-calculus model is process centric; there is one process for each component or observable 
product, leading to a much simpler model. Andrew Phillips shows the simplicity of the SPiM process 
model in contrast with the chemical reactions model EH . 

While ODEs are low level and close to chemical reactions, Petri Nets constitute a graphical represen- 
tation of chemical reactions. A SPiM model is a low level representation of individual protein behavior 
using send and receive communication primitives over a channel, and Kappa is a graphical high level 
abstraction of chemical reactions. 

4 High Level Notation 

Both ODEs and Petri Nets correspond closely to chemical reactions, and for the average biologist, they 
are relatively easy to understand. However, according to our experience in the classroom and in the lab, 
communicating with biology students and scientists, understanding the stochastic Pi-calculus, the Kappa 
language and Bio-PEPA (to a lesser extent) to be able to model molecular processes, requires a con- 
siderable initial effort. While Kappa has a user friendly graphical syntax available through Cellucidate. 
SPiM, on the other hand, still needs encoding in the stochastic Pi-calculus. To address this issue and 
enable modeling using descriptions directly obtained from biological processes, a narrative language has 
been proposed for SPiM lfT5l . Here we propose an alternative high level notation, and we show how to 
systematically translate it into SPiM programs and also chemical reactions. 

We write a, b, c, d, etc. for species, r for reaction rate constants, A for actions, and P for a process 
consisting of a possibly empty sequence of actions. Table [5] shows a list of actions representing familiar 
biological processes. bind(a, b, c, r) means two reagents a and b bind together to generate the prod- 
uct c at rate r. Similarly, dimerize(a, b, c, r) means that a and b form a dimer denoted by c at rate r. 
dissociate(a, b, c, r) is the opposite action to bind(a, b, c, r); a is broken down into b and c at rate r. 
activate(a, b, c, r) and activateAnddissociate(a, b, c, d, r) basically represent the same biological pro- 
cess where the reagent a becomes activated in the presence of another reagent b at rate r. The difference 
is that activate only generates the activated state product c while activateAnddissociate generates the 
activated product, and the product simultaneously dissociates into two products c and d. phosphory- 
late(a, b, c, r) represents phosphorylation at rate r, where a is the protein reagent, b is the kinase, and c 
is the product. hydrolyze(a, b, r) represents the hydrolysis at rate r, where a is the reagent and b is the 
product. degrade(a, r) means a degrades at rate r. Notice that bind(a, b, c, r) and dimerize(a, b, c, r) 
are basically the same operation, but in order to distinguish these two different biological processes, we 
have two names. In biology, dimerization is the binding which can generate a dimer. 

Table [5] also defines the translation function SPiM, that given an action in high level notation returns 
the corresponding SPiM process, and SPiM*, its homomorphic extension to processes. To aide read- 
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Figure 5: (a). Simulation result of ODEs model from Math Works, (b). Simulation result of process model from SPiM. (c). 
Simulation result of Petri Nets model from Cell Illustrator, (d). Simulation result of Kappa language from Cellucidate. (e). 
Simulation result of Bio-PEPA from Eclipse Plug-in. In all five cases, RL represents receptor bound to ligand, Ga represents 
the activated G-protein, and G represents the deactivated G-protein. The x axis is time in seconds, and the y axis is intensity of 
states. 



50 Computational Modeling for the Activation Cycle of G-proteins by G-protein-coupled Receptors 



Table 5: High Level Notation 



Action A 



Process P 



:= bind(a, b, c, r) 
activate(a, b, c, r) 
phosphorylate(a, b, c, r) 
degrade(a, r) 

:= A;P|() 



dimerize(a, b, c, r) 
activateAnddissociate(a, b, c, d, r) 
dissociate (a, b, c, r) 
hydrolyze(a, b, r) 



SPiM(bind(a, b, c, r)) = {a -+- ch@r c\b -^- ch@r 
SPiM(phosphorylate(a, b, c, r)) = {a W ch@r c\b ^ ch@r b} 
SPiM(dimerize(a, b, c, r)) = {a W ch@r 



c\b -^^ n ^' ()} = SPiM(dimerize(a, b, c, r)) where ch is new 

where ch is new 



\ch@r 



c\b -± lch@r ()} = SPiM(bind(a, b, c, r)) where ch is new 

c . b ^ch@r- 



SPiM(activate(a, b, c, r)) = {a 
SPiM(activateAnddissociate(a, b, c, d, r)) = {a 



b} 



Ach@r t 



d:b 



Jch@r 



SPiM(dissociate(a, b, c, r)) = {a ^ dda y @r 
SPiM(hydrolyze(a, b, r)= {a -^ dda y @r b} 
SPiM(degrade(a, r)= {a -^ dda y @r ()} 



b} 



b | c} 



where ch is new 
where ch is new 



SPiM* (0) = 

SPiM* (A; P) = SPiM(A) U SPiM* (P) 



Table 6: Translation of High Level Notation into Chemical Reactions 

Reaction (bind(a, b, c, r)) = a + b — ^ r c 

Reaction (phosphorylate(a, b, c, r))) = a + b^ r c + b 

Reaction (dimerize(a, b, c, r)) = a + b c 

Reaction (activate(a, b, c, r)) = a + b c + b 

Reaction(activateAnddissociate(a, b, c, d, r)) = a + b c + d + b 

Reaction(dissociate(a, b, c, r)) = a — ^ r b + c 

Reaction (hydrolyze(a, b, r))) = a — ^ r b 

Reaction(degrade(a, r)) = a null 

Reaction* (()) = 

Reaction* (A; P) = Reaction(A) U Reaction* (P) 



ability, instead of using stochastic Pi-calculus code, we use the equivalent graph representation of SPiM 
processes 0. Similarly, Table [6] defines the translation functions Reaction and Reaction* that translate 
processes in high level notation into chemical reactions. 

Going back to our G-Protein example in Fig. [TJ bind(Gd, Gbg, G, 1.0) corresponds to step 1 to 
2, bind(R, L, RL, 3.32e-6) to step 2 to 3, activate Anddissociate(G, RL, Ga, Gbg, 1.0e-5) to the 
activation in steps 2 to 3 and the dissociation in step 3 to 4, dissociate(RL, R, L, 0.01) represents step 
4 to 5 and hydrolyze(Ga, Gd, 0.11); degrade(R, 4e-4); degrade(RL, 4e-3) correspond to step 5 to 1 
completing the cycle. 

The translation of the high level notation of the G-protein cycle into SPiM is as follows: 
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SPiM*(bind(Gd, Gbg, G, 1.0); bind(R, L, RL, 3.32e-6); 

activate Anddissociate(G, RL, Ga, Gbg,1.0e-5); dissociate(RL, R, L, 0.01); 

hydrolyze(Ga, Gd, 0.11); degrade(R, 4e-4); degrade(RL, 4e-3)) 

{Gd ^-bindbm.o G . Gbg _^mndb@i.o ()} uSPiM*(bind(R, L, RL, 3.32e-6); 
activateAnddissociate(G, RL, Ga, Gbg,1.0e-5); dissociate(RL, R, L, 0.01); 
hydrolyze(Ga, Gd, 0.11); degrade(R, 4e-4); degrade(RL, 4e-3)) 



{ G d _^lbindb@1.0 q. Qfog _^Windb@l.O Q j j _^\bindR@ 3. 32e -6 .£ _^?bindR@3.32e-6 Q j j 
| G ^switch@1.0e-5 Gq | Q bg . RL ^\switch@1.0e-5 RL} U {RL ^ dela y @0M R | L}U 
{Ga -^delay@0.U GJ | y _^delay@4e-4 Q| y _^Z«y @4e-3 Q| 

Finally, the translation of the high level notation of the G-protein cycle into chemical reactions is as 
follows: 

Reaction*(bind(Gd, Gbg, G, 1.0); bind(R, L, RL, 3.32e-6); 
activateAnddissociate(G, RL, Ga, Gbg,1.0e-5); dissociate(RL, R, L, 0.01); 
hydrolyze(Ga, Gd, 0.11); degrade(R, 4e-4); degrade(RL, 4e-3)) 

Gd + Gbg -^ L0 GU Reaction*(bind(R, L, RL, 3.32e-6); 
activateAnddissociate(G, RL, Ga, Gbg,1.0e-5); dissociate(RL, R, L, 0.01); 
hydrolyze(Ga, Gd, 0.11); degrade(R, 4e-4); degrade(RL, 4e-3)) 



Gd + Gbg -> l G U R + L ^ 3 32 ^ 6 tfLU 

G + tfL^ 10 ^ 5 Ga + Gbg + RLURL^ 001 R + LU 

Ga -^°' n GdUR -^ 4e ~ 4 null URL -^ 4e ~ 3 null 

Kahramanogullari. et. al. fl31 , propose a narrative language that is translated into a model with three 
kinds of sentences (association, dissociation and transformation), where species have explicit binding 
sites (as in Kappa). They also show how to translate models into SPiM and how to obtain a model using 
the narrative language. Our high level notation is an alternative to their narrative language, however we 
translate it directly into SPiM, and we do not need a notion of model. In their example of Appendix A, 
the generated code for processes FcR2 and FcR3, for example, does not correspond to any species in the 
biological example being modeled. 

5 Conclusion 

Although building computational models of cellular processes is not new, to the best of our knowledge, 
this is the first time that a model has been implemented in five different formalisms. Furthermore, our 
classroom experience shows that this survey paper is a valuable tutorial introduction to bio-modeling 
in these diverse frameworks. In order to make a formal comparison for all five formalisms, one would 
need to define translation functions between the different formalisms and a notion of correctness for each 
translation. 
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Starting from an exiting ODEs model of the activation cycle of G-proteins by G-protein coupled 
receptors, we construct four simulations in four different formalisms: Stochastic Pi-Calculus, Stochastic 
Petri Nets, Kappa, and Bio-PEPA. We also show how to scale initial concentrations and reaction rates 
for this specific example to compensate for the fact that solving differential equations in MatLab is faster 
than executing a concurrent model. 

Finally, our high level notation is an abstraction for both SPiM and chemical reactions. Our high 
level notation corresponds to commonly-occurring reaction schemes easily identifiable by biologists, 
and it is an alternative to the narrative language proposed for SPiM. Our high level notation represents 
the G-protein cycle quite concisely, it is accessible to the average biologist, and it translates naturally 
into SPiM and chemical reactions, suggesting that our models could be effective instruments to assist in 
biomedical research. 

This notation naturally induces a typed modeling language that is currently being explored, where 
a,b h bind@r;P reduces to complex(a,b) h P, where complex is only defined for specific species, for 
example. 
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